Trait Ordination Plant Community Analysis

Trait Ordination Plant Community Analysis

1. traits

   1. above- and belowground

      1. height & SLA
         

   2. acquisitive vs. conservative resource-use strategies

      1. rooting depth

   3. collaborative vs. individual resource-use strategies

      1. fine root proportion

      2. thickness of roots

Goal: Examine how traits shift in community compositional (ordinational) space, placing each trait on an axis.

Author: Caryn D. Iwanaga

Updated: 02/06/2025

Libraries

library(tidyverse)
library(dplyr)
library(ggplot2)
library(httr) # read out Dropbox folders
library(vegan)
library(readxl)

Read in Data

cover.dat.labels <- read.csv("https://www.dropbox.com/scl/fi/up8nnzkcpchsr45f8cm92/Compost_Cover_LongClean.csv?rlkey=z2tvaj8t6khadef7ydz782zka&st=qwef9ys0&dl=1") %>%
  mutate(
    nut_trt = factor(ifelse(nut_trt =="C", "+compost", ifelse(nut_trt =="F", "+N fertilizer", ifelse(nut_trt =="N", "control", NA))), levels = c("control", "+N fertilizer", "+compost")),
    ppt_trt = ifelse(ppt_trt == "D", "dry", ifelse(ppt_trt == "XC", "control", ifelse(ppt_trt == "W", "wet", NA))))


# Define a reusable color scale
# custom_colors <- scale_fill_manual(
#   values = c(
#     "control" = rgb(195, 197, 193, maxColorValue = 250),
#     "+N fertilizer" = rgb(87, 62, 92, maxColorValue = 225),
#     "+compost" = rgb(167, 156, 109, maxColorValue = 225)
#   )
# )

custom_colors <- scale_fill_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

# Set a global theme
theme_set(theme_bw())

Make site x species matrix

1. Calculate relative abundances of species for each site
   1. confused how to aggregate abundances because even grouping by plotid (nut_trt and ppt_trt), year, species, and grazing history (block 1&2 – high, block 3&4 – low) and averaged, the sd are as high as 45.96!
      1. proceeded with these high sd to see what happens (ask Lauren)
   2. try how it is when you split everything by each block (more representative?)

site.df <- 
  cover.dat.labels[, c(1:7, 18:19)] %>% #subset columns
  mutate(
    grazing_hist = ifelse(block == 1 | block == 2, "high", ifelse(block == 3 | block == 4, "low", NA))
  ) %>% # separate by block (grazing intensities)
  group_by(nut_trt, ppt_trt, yr, code4, grazing_hist) %>% # make each code unique -- one value for species abundance for each trt
  summarise(
    # pct_cover,
    abundance = mean(pct_cover)
    # sd = sd(pct_cover)
  )

`summarise()` has grouped output by 'nut_trt', 'ppt_trt', 'yr', 'code4'. You can override using the `.groups` argument.

matrix0 <- spread(
  site.df, 
  code4, 
  abundance, 
  fill = 0)

# make rownames ----
rownames(matrix0) <- paste(matrix0$nut_trt, matrix0$ppt_trt, matrix0$yr, matrix0$grazing_hist, sep = "_")

Warning: Setting row names on a tibble is deprecated.

# delete columns ---- 
matrix <- matrix0[, c(5:79)]

# relativize data with Bray-Curtis dissimilarity
matrix.bray <- decostand(matrix, "total")
rownames(matrix.bray) <- rownames(matrix0)

Split matrix by each block as well

site.block.df <- 
  cover.dat.labels[, c(1:7, 18:19)] %>% #subset columns
  group_by(nut_trt, ppt_trt, yr, code4, block) %>% # make each code unique -- one value for species abundance for each trt
  summarise(
    # pct_cover,
    abundance = mean(pct_cover)
    # sd = sd(pct_cover)
  )

`summarise()` has grouped output by 'nut_trt', 'ppt_trt', 'yr', 'code4'. You can override using the `.groups` argument.

site.block.sd.df <- 
  cover.dat.labels[, c(1:7, 18:19)] %>% #subset columns
  group_by(nut_trt, ppt_trt, yr, code4, block) %>% # make each code unique -- one value for species abundance for each trt
  summarise(
    pct_cover,
    abundance = mean(pct_cover),
    sd = sd(pct_cover)
  )

Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in dplyr 1.1.0.
Please use `reframe()` instead.
When switching from `summarise()` to `reframe()`, remember that `reframe()` always returns an ungrouped data frame and adjust accordingly.`summarise()` has grouped output by 'nut_trt', 'ppt_trt', 'yr', 'code4', 'block'. You can override using the `.groups` argument.

matrix0.block <- spread(
  site.block.df, 
  code4, 
  abundance, 
  fill = 0)

# make rownames ----
rownames(matrix0.block) <- paste(matrix0.block$nut_trt, matrix0.block$ppt_trt, matrix0.block$yr, matrix0.block$block, sep = "_")

Warning: Setting row names on a tibble is deprecated.

# delete columns ---- 
matrix.block <- matrix0.block[, -c(0:4)]

# relativize data with Bray-Curtis dissimilarity
matrix.block.bray <- decostand(matrix.block, "total")
rownames(matrix.block.bray) <- rownames(matrix0.block)

Run NMDS

Rules of thumb:

• < 0.05 is excellent

• 0.05-0.1 is good

• 0.1-0.2 is fair

• 0.2-0.3 is cause for concern…

• > 0.3 is poor, random

All years – low & high grazing

Stress: fair

0.1966553 

Grouping by year and grazing history had most significant clustering shown.

metaMDS()

nmds.comp <- metaMDS(matrix.bray, k=2, trymax = 25)

Run 0 stress 0.2020736 
Run 1 stress 0.2044105 
Run 2 stress 0.2084765 
Run 3 stress 0.2084764 
Run 4 stress 0.2026778 
Run 5 stress 0.2018338 
... New best solution
... Procrustes: rmse 0.06193884  max resid 0.2323138 
Run 6 stress 0.1971149 
... New best solution
... Procrustes: rmse 0.02534655  max resid 0.1298511 
Run 7 stress 0.2027879 
Run 8 stress 0.20163 
Run 9 stress 0.2033608 
Run 10 stress 0.2017485 
Run 11 stress 0.2176904 
Run 12 stress 0.2013629 
Run 13 stress 0.1966668 
... New best solution
... Procrustes: rmse 0.01514645  max resid 0.1046137 
Run 14 stress 0.1991103 
Run 15 stress 0.202025 
Run 16 stress 0.2033208 
Run 17 stress 0.2020741 
Run 18 stress 0.2034839 
Run 19 stress 0.2043941 
Run 20 stress 0.1991947 
Run 21 stress 0.2128964 
Run 22 stress 0.1991946 
Run 23 stress 0.2039391 
Run 24 stress 0.2133823 
Run 25 stress 0.196656 
... New best solution
... Procrustes: rmse 0.004511995  max resid 0.02760931 
*** Best solution was not repeated -- monoMDS stopping criteria:
    24: stress ratio > sratmax
     1: scale factor of the gradient < sfgrmin

nmds.comp # pretty not great stress

Call:
metaMDS(comm = matrix.bray, k = 2, trymax = 25) 

global Multidimensional Scaling using monoMDS

Data:     matrix.bray 
Distance: bray 

Dimensions: 2 
Stress:     0.196656 
Stress type 1, weak ties
Best solution was not repeated after 25 tries
The best solution was from try 25 (random start)
Scaling: centring, PC rotation, halfchange scaling 
Species: expanded scores based on ‘matrix.bray’ 

stressplot(nmds.comp)

plot(nmds.comp, type="t")

nmds1 <- as.data.frame(scores(nmds.comp, choices=c(1), display=c("sites"))) %>%
  mutate(matrix0[, c(0:4)])
nmds2 <- as.data.frame(scores(nmds.comp, choices=c(2), display=c("sites")))

nmds_dat <- cbind(nmds1, nmds2) %>%
  select(nut_trt:grazing_hist, NMDS1, NMDS2)

plot

ggplot(nmds_dat, aes(NMDS1, NMDS2,  fill = nut_trt, color = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Nutrient Treatment") +
  custom_colors +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
Please use `linewidth` instead.

ggplot(nmds_dat, aes(NMDS1, NMDS2, color = as.factor(yr), fill = as.factor(yr))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Year")

ggplot(nmds_dat, aes(NMDS1, NMDS2, fill = as.factor(ppt_trt), color = as.factor(ppt_trt))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Precipitation Treatment") +
  custom_colors +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

ggplot(nmds_dat, aes(NMDS1, NMDS2,  fill = as.factor(grazing_hist), color = as.factor(grazing_hist))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Grazing History")

ggplot(nmds_dat, aes(NMDS1, NMDS2,  color = as.factor(grazing_hist), fill = as.factor(yr))) +
  geom_point(aes(fill = as.factor(yr)), shape = 21, size = 4, stroke = 1) +
  stat_ellipse(aes(color = as.factor(yr)), 
                geom = "polygon",
               fill = NA,
               size = 1) +
  labs(title = "Grouped by Grazing History & Year") +
  scale_color_manual(
  values = c(
    "low" = "cyan",
    "high" = "brown",
  "2019" = "salmon",
  "2020" = "limegreen",
  "2021" = "cornflowerblue")
  )

  
ggplot(nmds_dat, aes(NMDS1, NMDS2, color = nut_trt, fill = as.factor(yr))) +
  geom_point(shape = 21, size = 4, stroke = 1) +
  stat_ellipse(aes(color = as.factor(nut_trt)),
             geom = "polygon",
             fill = NA,  # Makes the ellipse transparent inside
             size = 1) +
  labs(title = "Grouped by Nutrient Treatment & Year") +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

All years – blocks

Stress: cause for concern

0.2329306

metaMDS()

nmds.block.comp <- metaMDS(matrix.block.bray, k=2, trymax = 25)
nmds.block.comp
stressplot(nmds.block.comp)
plot(nmds.block.comp, type="t")

nmds1 <- as.data.frame(scores(nmds.block.comp, choices=c(1), display=c("sites"))) %>%
  mutate(matrix0.block[, c(0:4)])
nmds2 <- as.data.frame(scores(nmds.block.comp, choices=c(2), display=c("sites")))

nmds_dat.block <- cbind(nmds1, nmds2) %>%
  select(nut_trt:block, NMDS1, NMDS2)

plot

ggplot(nmds_dat.block, aes(NMDS1, NMDS2,  fill = nut_trt, color = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Nutrient Treatment - Blocks") +
  custom_colors +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

ggplot(nmds_dat.block, aes(NMDS1, NMDS2, color = as.factor(yr), fill = as.factor(yr))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Year - Blocks")

ggplot(nmds_dat.block, aes(NMDS1, NMDS2, fill = as.factor(ppt_trt), color = as.factor(ppt_trt))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Precipitation Treatment - Blocks") +
  custom_colors +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

ggplot(nmds_dat.block, aes(NMDS1, NMDS2,  fill = as.factor(block), color = as.factor(block))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Block")

ggplot(nmds_dat.block, aes(NMDS1, NMDS2,  color = as.factor(block), fill = as.factor(yr))) +
  geom_point(aes(fill = as.factor(yr)), shape = 21, size = 4, stroke = 1) +
  stat_ellipse(aes(color = as.factor(block)), 
                geom = "polygon",
               fill = NA,
               size = 1) +
  labs(title = "Grouped by Block & Year - Blocks") +
  scale_color_manual(
  values = c(
    # "low" = "cyan",
    # "high" = "brown",
  "2019" = "salmon",
  "2020" = "limegreen",
  "2021" = "cornflowerblue")
  )
  
ggplot(nmds_dat.block, aes(NMDS1, NMDS2, color = nut_trt, fill = as.factor(yr))) +
  geom_point(shape = 21, size = 4, stroke = 1) +
  stat_ellipse(aes(color = as.factor(nut_trt)),
             geom = "polygon",
             fill = NA,  # Makes the ellipse transparent inside
             size = 1) +
  labs(title = "Grouped by Nutrient Treatment & Year - Blocks") +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

All years - LOW GRAZING ONLY

Stress: fair

0.1869068

data wrangling

# relativize data with Bray-Curtis dissimilarity
matrix.low <- matrix0 %>%
  filter(grazing_hist == "low")

row.low <- paste(matrix.low$nut_trt, matrix.low$ppt_trt, matrix.low$yr, matrix.low$grazing_hist, sep = "_")

matrix.low <- matrix.low[, c(5:79)]

rownames(matrix.low) <- row.low

matrix.bray.low <- decostand(matrix.low, "total")

metaMDS()

nmds.comp.low <- metaMDS(matrix.bray.low, k=2, trymax = 25)
nmds.comp.low # pretty not great stress
stressplot(nmds.comp.low)
plot(nmds.comp.low, type="t")

nmds1 <- as.data.frame(scores(nmds.comp.low, choices=c(1), display=c("sites"))) %>%
  mutate(
    nut_trt = sapply(strsplit(rownames(matrix.low), "_"), `[`, 1),
    ppt_trt = sapply(strsplit(rownames(matrix.low), "_"), `[`, 2),
    yr = sapply(strsplit(rownames(matrix.low), "_"), `[`, 3),
    grazing_hist = sapply(strsplit(rownames(matrix.low), "_"), `[`, 4)
  )
nmds2 <- as.data.frame(scores(nmds.comp.low, choices=c(2), display=c("sites")))

nmds_low_dat <- cbind(nmds1, nmds2) %>%
  select(nut_trt:grazing_hist, NMDS1, NMDS2)

plot

ggplot(nmds_low_dat, aes(NMDS1, NMDS2,  color = as.factor(yr), shape = nut_trt)) +
  geom_point(size = 3, stroke = 1) + 
  labs(title = "2019-2021 (LOW) - Grouping by Year & Nutrient Treatment")

ggplot(nmds_low_dat, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  labs(title = "2019-2021 (LOW) - Grouping by Nutrient Treatment")+ 
  custom_colors


ggplot(nmds_low_dat, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2019-2021 (LOW) - Grouping by Precipitation Treatment") 

# add year ellipses ---- 
ggplot(nmds_low_dat, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(aes(color = as.factor(yr)), 
                geom = "polygon",
               fill = NA,
               size = 1) +
  labs(title = "2019-2021 (LOW) - Grouping by Nutrient Treatment")+ 
  custom_colors


ggplot(nmds_low_dat, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
    stat_ellipse(aes(color = as.factor(yr)), 
                geom = "polygon",
               fill = NA,
               size = 1) +
  custom_colors + 
  labs(title = "2019-2021 (LOW) - Grouping by Precipitation Treatment") 

Make each individual year NMDS to see if there’s better stress

2019

data wrangling

# relativize data with Bray-Curtis dissimilarity
matrix.2019 <- matrix0 %>%
  filter(yr == 2019)

row.2019 <- paste(matrix.2019$nut_trt, matrix.2019$ppt_trt, matrix.2019$yr, matrix.2019$grazing_hist, sep = "_")

matrix.2019 <- matrix.2019[, c(5:79)]

rownames(matrix.2019) <- row.2019

matrix.bray.2019 <- decostand(matrix.2019, "total")

nmds 2019

Stress: Excellent

0.0560178

metaMDS()

nmds.2019 <- metaMDS(matrix.bray.2019, k=2, trymax = 25)
nmds.2019 # pretty not great stress
stressplot(nmds.2019)
plot(nmds.2019, type="t")

nmds1.2019 <- as.data.frame(scores(nmds.2019, choices=c(1), display=c("sites"))) %>%
  mutate(
    nut_trt = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 1),
    ppt_trt = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 2),
    yr = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 3),
    grazing_hist = sapply(strsplit(rownames(matrix.2019), "_"), `[`, 4)
  )
nmds2.2019 <- as.data.frame(scores(nmds.2019, choices=c(2), display=c("sites")))

nmds_2019_dat <- cbind(nmds1.2019, nmds2.2019) %>%
  select(nut_trt:grazing_hist, NMDS1, NMDS2)

plot

ggplot(nmds_2019_dat, aes(NMDS1, NMDS2,  color = grazing_hist, shape = nut_trt)) +
  geom_point(size = 3, stroke = 1) + 
  labs(title = "2019 - Grouping by Grazing History & Nutrient Treatment")

ggplot(nmds_2019_dat, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  labs(title = "2019 - Grouping by Nutrient Treatment")+ 
  custom_colors


ggplot(nmds_2019_dat, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2019 - Grouping by Precipitation Treatment") 

2020

data wrangling

# relativize data with Bray-Curtis dissimilarity
matrix.2020 <- matrix0 %>%
  filter(yr == 2020)

row.2020 <- paste(matrix.2020$nut_trt, matrix.2020$ppt_trt, matrix.2020$yr, matrix.2020$grazing_hist, sep = "_")

matrix.2020 <- matrix.2020[, c(5:79)]

rownames(matrix.2020) <- row.2020

matrix.bray.2020 <- decostand(matrix.2020, "total")

nmds 2020

Stress: fair

0.1522495

metaMDS()

nmds.2020 <- metaMDS(matrix.bray.2020, k=2, trymax = 25)
nmds.2020
stressplot(nmds.2020)
plot(nmds.2020, type="t")

nmds1.2020 <- as.data.frame(scores(nmds.2020, choices=c(1), display=c("sites"))) %>%
  mutate(
    nut_trt = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 1),
    ppt_trt = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 2),
    yr = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 3),
    grazing_hist = sapply(strsplit(rownames(matrix.2020), "_"), `[`, 4)
  )
nmds2.2020 <- as.data.frame(scores(nmds.2020, choices=c(2), display=c("sites")))

nmds_2020_dat <- cbind(nmds1.2020, nmds2.2020) %>%
  select(nut_trt:grazing_hist, NMDS1, NMDS2)

plot

ggplot(nmds_2020_dat, aes(NMDS1, NMDS2,  color = grazing_hist, shape = nut_trt)) +
  geom_point(size = 4, stroke = 1) + 
  labs(title = "2020 - Grouping by Grazing History & Nutrient Treatment") 

ggplot(nmds_2020_dat, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2020 - Grouping by Nutrient Treatment")

ggplot(nmds_2020_dat, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2020 - Grouping by Precipitation Treatment")

2021

data wrangling

# relativize data with Bray-Curtis dissimilarity
matrix.2021 <- matrix0 %>%
  filter(yr == 2021)

row.2021 <- paste(matrix.2021$nut_trt, matrix.2021$ppt_trt, matrix.2021$yr, matrix.2021$grazing_hist, sep = "_")

matrix.2021 <- matrix.2021[, c(5:79)]

rownames(matrix.2021) <- row.2021

matrix.bray.2021 <- decostand(matrix.2021, "total")

2021 - LOW

Stress

0.08067789

# relativize data with Bray-Curtis dissimilarity
matrix.2021.1 <- matrix0 %>%
  filter(yr == 2021,
         grazing_hist == "low") 

row.2021.1 <- paste(matrix.2021.1$nut_trt, matrix.2021.1$ppt_trt, matrix.2021.1$yr, matrix.2021.1$grazing_hist, sep = "_")

matrix.2021.1 <- matrix.2021.1[, c(5:79)]

rownames(matrix.2021.1) <- row.2021.1

matrix.bray.2021.1 <- decostand(matrix.2021.1, "total")

nmds.2021.1 <- metaMDS(matrix.bray.2021.1, k=2, trymax = 25)
nmds.2021.1
stressplot(nmds.2021.1)
plot(nmds.2021.1, type="t")
dummy <- as.data.frame(scores(nmds.2021.1, choices=c(1), display=c("sites")))

nmds1.2021.1.df <- as.data.frame(scores(nmds.2021.1, choices=c(1), display=c("sites"))) %>%
  mutate(trt = rownames(dummy)) %>%
  separate(trt, c("nut_trt", "ppt_trt", "pr", "grazing_hist"), sep="_")

nmds2.2021.1 <- as.data.frame(scores(nmds.2021.1, choices=c(2), display=c("sites")))

nmds_2021_dat.1 <- cbind(nmds1.2021.1.df, nmds2.2021.1) %>%
  select(nut_trt:grazing_hist, NMDS1, NMDS2)

# ggplot(nmds_2021_dat.1, aes(NMDS1, NMDS2,  color = grazing_hist, shape = nut_trt)) +
#   geom_point(size = 4, stroke = 0.5) +
#   labs(title = "2021 - Grouping by Grazing History & Nutrient Treatment")

ggplot(nmds_2021_dat.1, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (LOW) - Grouping by Nutrient Treatment")

ggplot(nmds_2021_dat.1, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (LOW) - Grouping by Precipitation Treatment")

2021 - HIGH

Stress

0.08065034

# relativize data with Bray-Curtis dissimilarity
matrix.2021.2 <- matrix0 %>%
  filter(yr == 2021,
         grazing_hist == "high") 

row.2021.2 <- paste(matrix.2021.2$nut_trt, matrix.2021.2$ppt_trt, matrix.2021.2$yr, matrix.2021.2$grazing_hist, sep = "_")

matrix.2021.2 <- matrix.2021.2[, c(5:79)]

rownames(matrix.2021.2) <- row.2021.2

matrix.bray.2021.2 <- decostand(matrix.2021.2, "total")

nmds.2021.2 <- metaMDS(matrix.bray.2021.2, k=2, trymax = 25)
nmds.2021.2
stressplot(nmds.2021.2)
plot(nmds.2021.2, type="t")
dummy <- as.data.frame(scores(nmds.2021.2, choices=c(1), display=c("sites")))

nmds1.2021.2.df <- as.data.frame(scores(nmds.2021.2, choices=c(1), display=c("sites"))) %>%
  mutate(trt = rownames(dummy)) %>%
  separate(trt, c("nut_trt", "ppt_trt", "pr", "grazing_hist"), sep="_")

nmds2.2021.2 <- as.data.frame(scores(nmds.2021.2, choices=c(2), display=c("sites")))

nmds_2021_dat.2 <- cbind(nmds1.2021.2.df, nmds2.2021.2) %>%
  select(nut_trt:grazing_hist, NMDS1, NMDS2)

ggplot(nmds_2021_dat.2, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (HIGH) - Grouping by Nutrient Treatment")

ggplot(nmds_2021_dat.2, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (HIGH) - Grouping by Precipitation Treatment")

nmds 2021

Stress: fair

0.1230831

metaMDS()

nmds.2021 <- metaMDS(matrix.bray.2021, k=2, trymax = 25)
nmds.2021
stressplot(nmds.2021)
plot(nmds.2021, type="t")

nmds1.2021 <- as.data.frame(scores(nmds.2021, choices=c(1), display=c("sites"))) %>%
  mutate(
    nut_trt = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 1),
    ppt_trt = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 2),
    yr = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 3),
    grazing_hist = sapply(strsplit(rownames(matrix.2021), "_"), `[`, 4)
  )
nmds2.2021 <- as.data.frame(scores(nmds.2021, choices=c(2), display=c("sites")))

nmds_2021_dat <- cbind(nmds1.2021, nmds2.2021) %>%
  select(nut_trt:grazing_hist, NMDS1, NMDS2)

plot

ggplot(nmds_2021_dat, aes(NMDS1, NMDS2,  color = grazing_hist, shape = nut_trt)) +
  geom_point(size = 4, stroke = 0.5) +
  labs(title = "2021 - Grouping by Grazing History & Nutrient Treatment")

ggplot(nmds_2021_dat, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 - Grouping by Nutrient Treatment")

ggplot(nmds_2021_dat, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 - Grouping by Precipitation Treatment")

2021 — blocks

data wrangling

# relativize data with Bray-Curtis dissimilarity
matrix.2021.block <- matrix0.block %>%
  filter(yr == 2021)

row.2021.block <- paste(matrix.2021.block$nut_trt, matrix.2021.block$ppt_trt, matrix.2021.block$yr, matrix.2021.block$block, sep = "_")

matrix.2021.block <- matrix.2021.block[, c(5:79)]

rownames(matrix.2021.block) <- row.2021.block

matrix.bray.block.2021 <- decostand(matrix.2021.block, "total")

nmds 2021

Stress: fair

0.1702944

metaMDS()

nmds.2021.block <- metaMDS(matrix.bray.block.2021, k=2, trymax = 25)
nmds.2021.block
stressplot(nmds.2021.block)
plot(nmds.2021.block, type="t")

nmds1.2021.block <- as.data.frame(scores(nmds.2021.block, choices=c(1), display=c("sites"))) %>%
  mutate(
    nut_trt = sapply(strsplit(rownames(matrix.2021.block), "_"), `[`, 1),
    ppt_trt = sapply(strsplit(rownames(matrix.2021.block), "_"), `[`, 2),
    yr = sapply(strsplit(rownames(matrix.2021.block), "_"), `[`, 3),
    block = sapply(strsplit(rownames(matrix.2021.block), "_"), `[`, 4)
  )
nmds2.2021.block <- as.data.frame(scores(nmds.2021.block, choices=c(2), display=c("sites")))

nmds_2021_dat.block <- cbind(nmds1.2021.block, nmds2.2021.block) %>%
  select(nut_trt:block, NMDS1, NMDS2)

plot

ggplot(nmds_2021_dat.block, aes(NMDS1, NMDS2,  color = block, shape = nut_trt)) +
  geom_point(size = 4, stroke = 0.5) +
  labs(title = "2021 (Block) - Grouping by Block & Nutrient Treatment")

ggplot(nmds_2021_dat.block, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (Block) - Grouping by Nutrient Treatment")

ggplot(nmds_2021_dat.block, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (Block) - Grouping by Precipitation Treatment")

2021 — blocks (LOW)

data wrangling

# relativize data with Bray-Curtis dissimilarity
matrix.2021.block.low <- matrix0.block %>%
  filter(yr == 2021,
         block == 3 | block == 4)

row.2021.block.low <- paste(matrix.2021.block.low$nut_trt, matrix.2021.block.low$ppt_trt, matrix.2021.block.low$yr, matrix.2021.block.low$block, sep = "_")

matrix.2021.block.low <- matrix.2021.block.low[, -c(0:4)]

rownames(matrix.2021.block.low) <- row.2021.block.low

matrix.bray.block.2021.low <- decostand(matrix.2021.block.low, "total")

nmds 2021

Stress: fair

0.1282296

metaMDS()

nmds.2021.block.low <- metaMDS(matrix.bray.block.2021.low, k=2, trymax = 25)
nmds.2021.block.low
stressplot(nmds.2021.block.low)
plot(nmds.2021.block.low, type="t")

nmds1.2021.block.low <- as.data.frame(scores(nmds.2021.block.low, choices=c(1), display=c("sites"))) %>%
  mutate(
    nut_trt = sapply(strsplit(rownames(matrix.2021.block.low), "_"), `[`, 1),
    ppt_trt = sapply(strsplit(rownames(matrix.2021.block.low), "_"), `[`, 2),
    yr = sapply(strsplit(rownames(matrix.2021.block.low), "_"), `[`, 3),
    block = sapply(strsplit(rownames(matrix.2021.block.low), "_"), `[`, 4)
  )
nmds2.2021.block.low <- as.data.frame(scores(nmds.2021.block.low, choices=c(2), display=c("sites")))

nmds_2021_dat.block.low <- cbind(nmds1.2021.block.low, nmds2.2021.block.low) %>%
  select(nut_trt:block, NMDS1, NMDS2)

plot

ggplot(nmds_2021_dat.block.low, aes(NMDS1, NMDS2,  color = block, shape = nut_trt)) +
  geom_point(size = 4, stroke = 0.5) +
  labs(title = "2021 (Block; LOW) - Grouping by Block & Nutrient Treatment")

ggplot(nmds_2021_dat.block.low, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (Block; LOW) - Grouping by Nutrient Treatment")

ggplot(nmds_2021_dat.block.low, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (Block; LOW) - Grouping by Precipitation Treatment")

2021 — blocks (HIGH)

data wrangling

# relativize data with Bray-Curtis dissimilarity
matrix.2021.block.high <- matrix0.block %>%
  filter(yr == 2021,
         block == 1 | block == 2)

row.2021.block.high <- paste(matrix.2021.block.high$nut_trt, matrix.2021.block.high$ppt_trt, matrix.2021.block.high$yr, matrix.2021.block.high$block, sep = "_")

matrix.2021.block.high <- matrix.2021.block.high[, -c(0:4)]

rownames(matrix.2021.block.high) <- row.2021.block.high

matrix.bray.block.2021.high <- decostand(matrix.2021.block.high, "total")

nmds 2021

Stress: fair

0.1412997

metaMDS()

nmds.2021.block.high <- metaMDS(matrix.bray.block.2021.high, k=2, trymax = 25)
nmds.2021.block.high
stressplot(nmds.2021.block.high)
plot(nmds.2021.block.high, type="t")

nmds1.2021.block.high <- as.data.frame(scores(nmds.2021.block.high, choices=c(1), display=c("sites"))) %>%
  mutate(
    nut_trt = sapply(strsplit(rownames(matrix.2021.block.high), "_"), `[`, 1),
    ppt_trt = sapply(strsplit(rownames(matrix.2021.block.high), "_"), `[`, 2),
    yr = sapply(strsplit(rownames(matrix.2021.block.high), "_"), `[`, 3),
    block = sapply(strsplit(rownames(matrix.2021.block.high), "_"), `[`, 4)
  )
nmds2.2021.block.high <- as.data.frame(scores(nmds.2021.block.high, choices=c(2), display=c("sites")))

nmds_2021_dat.block.high <- cbind(nmds1.2021.block.high, nmds2.2021.block.high) %>%
  select(nut_trt:block, NMDS1, NMDS2)

plot

ggplot(nmds_2021_dat.block.high, aes(NMDS1, NMDS2,  color = block, shape = nut_trt)) +
  geom_point(size = 4, stroke = 0.5) +
  labs(title = "2021 (Block; HIGH) - Grouping by Block & Nutrient Treatment")

ggplot(nmds_2021_dat.block.high, aes(NMDS1, NMDS2,  fill = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (Block; HIGH) - Grouping by Nutrient Treatment")

ggplot(nmds_2021_dat.block.high, aes(NMDS1, NMDS2,  fill = ppt_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  custom_colors + 
  labs(title = "2021 (Block; HIGH) - Grouping by Precipitation Treatment")

Make species x trait matrix

Greenhouse traits

• height, leaf dry matter content, specific leaf area, root mass fraction, relative growth rate, coarse root diameter, root density, specific root length of fine roots, specific root length of coarse roots

Field traits

• height, leaf dry matter content, specific leaf area, seed mass, leaf C, leaf N

“H1: Pairing field results with a greenhouse 15 N experiment will allow us to relate species responses to their ability to uptake organic N”

“In addition, we will conduct a greenhouse experiment to isolate the mechanism (mineralized N vs organic N uptake rates ) by which compost may affect species composition”

Workflow plan:

• overlay species trait data as vectors on the NMDS

  • only greenhouse traits available (only sla data from field)

    • use PC scores as composite trait axes OR

    • average the traits for each species

Read in trait data

rda(traits.gh, scale = TRUE)

Error in colMeans(x, na.rm = TRUE) : 'x' must be numeric

Investigate trait data, looking in Greenhouse trait folder

Trait data is identical. Will clean data in

# dropbox.main <- read.csv("https://www.dropbox.com/scl/fi/c8vk31a807xg48jez1j03/GreenhouseTraits.csv?rlkey=ca736ea27k79mo7h8dbn67qiq&st=0h1bs53j&dl=1")
# 
# dropbox.gh <- read.csv("https://www.dropbox.com/scl/fi/67aaubdydiy4ixffmwvfe/GreenhouseTraits.csv?rlkey=rm2eqeqome7hgoedhnal0gx6d&st=vkym6evx&dl=1")
# 
# all.equal(dropbox.main, dropbox.gh)
# identical(dropbox.main, dropbox.gh)

Summarize trait data (GH) for each species

means

# make trait species ID match matrix data (notes shown below) -------
traits.cleaned <- traits.cleaned %>%
  mutate(
    ID = ifelse(
      ID == "Agoseris", "AGSP", ifelse(
        ID == "AVEBAR", "AVBA", ifelse(
          ID == "BRDR", "BRDI", ID
        )
      )
    )
  )

# summary(traits.cleaned)
traits.means <- traits.cleaned %>%
  group_by(ID) %>%
  summarise(
    across(4:(ncol(traits.gh) - 1), list(mean = ~mean(.x, na.rm = TRUE)
  )))
  

# make site by species matrix match traits species
matrix.names <- data.frame(colnames(matrix.bray)) #species comp (75 count)
trait.names <- data.frame(traits.means$ID) #trait species (80 count)

intersect.names <- intersect(colnames(matrix.bray), traits.means$ID) # same species (43 count)
setdiff(colnames(matrix.bray), traits.means$ID)

 [1] "ASSP"  "AST1"  "AST2"  "AST3"  "CIQU"  "CRTI"  "CYEC"  "DACA"  "DIMU"  "GAGR"  "GAPA"  "GAPH"  "HYSP"  "JUBU"  "LYHY" 
[16] "MAD1"  "MAD2"  "NAPU"  "PAVI"  "SABI1" "SABI2" "SIGA"  "THGR"  "TRSP"  "UNBU"  "UNF1"  "UNF3"  "UNF4"  "UNF8"  "UNFR" 
[31] "UNGR"  "XAST" 

setdiff(traits.means$ID, colnames(matrix.bray))

 [1] "ACMI"      "ANARp"     "AVFA"      "BRCA"      "BRHOp"     "BRNIf"     "BRNIp"     "CLPUf"     "CLPUp"     "CYDA"     
[11] "Crassula"  "ESCA"      "FEMI"      "FEMY"      "GITRf"     "GITRp"     "HYGL"      "LACA"      "LENIf"     "LENIp"    
[21] "LOPU"      "LUBI"      "Leontodon" "MAELf"     "MAELp"     "MICA"      "NEMA"      "PLERf"     "PLERp"     "RUPU"     
[31] "Sanicula"  "TRCI"      "TRHIpi"    "TRWIf"     "TRWIp"     "TRWIpi"    "VIVA"     

# check matrix codes to see if same species are named different things
unique.matrix.species <- cover.dat %>%
  dplyr::select(code4, species) %>%
  summarise(
    code4 = unique(code4),
    species = unique(species)
  )

Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in dplyr 1.1.0.
Please use `reframe()` instead.
When switching from `summarise()` to `reframe()`, remember that `reframe()` always returns an ungrouped data frame and adjust accordingly.

unique.trait.species <- traits.cleaned %>%
  dplyr::select(Taxon, ID) %>%
  distinct(ID, .keep_all = TRUE)
  

• change Trait name: ‘Agoseris’ –> ‘Agoseris unknown’; assuming the same as ‘AGSP’ –> ‘Agoseris sp.’ in matrix data (change to matrix name)

• What’s the difference between:

  • ‘ANAR’ and ‘ANARp’ in trait data?

  • ‘BRHO’ and ‘BRHOp’ in trait data?

  • ‘BRNIf’ and ‘BRNIp’ in trait data? (not in matrix data)

  • ‘CLPUf’ and ‘CLPUp’ in trait data? (not in matrix data)

  • ‘GITRf’ and ‘GITRp’ in trait data? (not in matrix data)

  • ‘LENIf’ and ‘LENIp’ in trait data? (not in matrix data)

  • ‘TRHI’ and ‘TRHIpi’ (same sp name) in trait data? (TRHI in matrix data)

  • ‘TRWIf’ and ‘TRWIp’ and ‘TRWIpi’ in trait data? (not in matrix data)

  • ‘Crassula’ and ‘Crassula tillaea’ (could they be the same?)

  • ‘Sanicula’ (unknown) and ‘SABI1’ or ‘SABI2’ (could they be the same?)

  • ‘MAELf’ (trait data) and ‘MAD1’ (Madia sp. 1 (“tall tarweed”)’ (could they be the same?)

  • ‘MAELp’ (trait data) and ‘MAD2’ (Madia sp. 1 (“small tarweed”)’ (could they be the same?)

  • ‘Trifolium eriocephalum’ (TRER) and ‘Triphysaria eriantha’ (TRER) (could they be the same?)

• change ‘AVEBAR’ (trait data) to ‘AVBA’ (matrix data)

• change BRDR –> BRDI

Overlay Trait Vectors onto NMDS

ALL YEARS – LOW AND HIGH GRAZING

Prepare NMDS and Traits

Stress:

0.2011127

# trait vectors onto the NMDS
envfit(nmds.comp.trait, traits.means.comp, permutations = 999)

Warning: longer object length is not a multiple of shorter object lengthError in vectorfit(X, P, permutations, strata, choices, w = w, ...) : 
  input data have non-matching numbers of observations

Plot

ggplot(nmds_dat_trait, aes(NMDS1, NMDS2,  fill = nut_trt, color = nut_trt)) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Nutrient Treatment") +
  custom_colors +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, color = as.factor(yr), fill = as.factor(yr))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Year")

ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, fill = as.factor(ppt_trt), color = as.factor(ppt_trt))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Precipitation Treatment") +
  custom_colors +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

ggplot(nmds_dat_trait, aes(NMDS1, NMDS2,  fill = as.factor(grazing_hist), color = as.factor(grazing_hist))) +
  geom_point(shape = 21, size = 4, stroke = 0.5) + 
  stat_ellipse(type = "norm", size = 1, alpha = 0.5) + # Ellipses for each group
  labs(title = "Grouped by Grazing History")

ggplot(nmds_dat_trait, aes(NMDS1, NMDS2,  color = as.factor(grazing_hist), fill = as.factor(yr))) +
  geom_point(aes(fill = as.factor(yr)), shape = 21, size = 4, stroke = 1) +
  stat_ellipse(aes(color = as.factor(yr)), 
                geom = "polygon",
               fill = NA,
               size = 1) +
  labs(title = "Grouped by Grazing History & Year") +
  scale_color_manual(
  values = c(
    "low" = "cyan",
    "high" = "brown",
  "2019" = "salmon",
  "2020" = "limegreen",
  "2021" = "cornflowerblue")
  )

  
ggplot(nmds_dat_trait, aes(NMDS1, NMDS2, color = nut_trt, fill = as.factor(yr))) +
  geom_point(shape = 21, size = 4, stroke = 1) +
  stat_ellipse(aes(color = as.factor(nut_trt)),
             geom = "polygon",
             fill = NA,  # Makes the ellipse transparent inside
             size = 1) +
  labs(title = "Grouped by Nutrient Treatment & Year") +
  scale_color_manual(
  values = c(
    "control" = rgb(220, 220, 220, maxColorValue = 255),      # Lighter Gray
    "+N fertilizer" = rgb(140, 40, 150, maxColorValue = 255), # Vibrant Purple
    "+compost" = rgb(220, 180, 60, maxColorValue = 255),       # Bright Golden Yellow
    "wet" = rgb(0, 0, 255, maxColorValue = 255),
    "dry" = rgb(183, 65, 14, maxColorValue = 255)
  )
)

# setdiff(colnames(matrix.bray), traits.gh.means$ID)
# setdiff(traits.gh.means$ID, colnames(matrix.bray))
traits.gh.means2 <- traits.gh.means %>%
  filter(ID %in% intersect(colnames(matrix.bray), traits.gh.means$ID))
rownames(traits.gh.means2) <- as.character(traits.gh.means2$ID)

matrix.bray %>%
  select(all_of(intersect(colnames(matrix.bray), traits.gh.means$ID))) # is this appropriate? should i relativize it all only with these species instead of just cutting out species?
nmds.comp.2 <- metaMDS(matrix.bray %>%
  select(all_of(intersect(colnames(matrix.bray), traits.gh.means$ID))), k=2, trymax = 25)
nmds.comp.2 # pretty not great stress
stressplot(nmds.comp.2)
plot(nmds.comp.2, type="t")

envfit(nmds.comp.2, traits.gh.means2, permutations = 999)